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ABSTRACT 

We develop an algorithm of separating the E and B modes of the CMB 
polarization from the noisy and discretized maps of Stokes parameter Q and 
t/ in a finite area. A key step of the algorithm is to take a wavelet- Galerkin 
discretization of the differential relation between the E, B and Q, U fields. This 
discretization allows derivative operator to be represented by a matrix, which is 
exactly diagonal in scale space, and narrowly banded in spatial space. We show 
that the effect of boundary can be eliminated by dropping a few DWT modes 
located on or nearby the boundary. This method reveals that the derivative 
operators will cause large errors in the E and B power spectra on small scales 
if the Q and U maps contain Gaussian noise. It also reveals that if the Q and 
U maps are random, these fields lead to the mixing of the E and B modes. 
Consequently, the B mode will be contaminated if the powers of E modes are 
much larger than that of B modes. Nevertheless, numerical tests show that the 
power spectra of both E and B on scales larger than the finest scale by a factor 
of 4 and higher can reasonably be recovered, even when the power ratio of E- to 
B-modes is as large as about 10^, and the signal-to- noise ratio is equal to 10 and 
higher. This is because the Galerkin discretization is free of false correlations, and 
keeps the contamination under control. As wavelet variables contain information 
of both spatial and scale spaces, the developed method is also effective to recover 
the spatial structures of the E and B mode fields. 
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1. Introduction 

The scalar component of primordial perturbations of the universe can be detected by 
the maps of temperature fluctuations of Cosmic Microwave Background Radiation (CMBR), 
while the tensor component of the primordial perturbations has to be probed by the maps of 
the Stokes parameter Q and U of the hnear polarization of the CMBR. A tensor field generally 
contains electric-like £^-mode and magnetic-like S-modes. In the linear regime, vortical mode 
of primordial perturbations do not grow during the clustering of density field, and therefore, 
the perturbed field initially has to be curl-free. That is, the primordial perturbations can 
only yield the £^-mode, but not 5-mode of the CMBR polarization field. On the other 
hand, B-mode perturbations can be produced by gravitational waves. Therefore, extracting 
the 5-mode information from CMBR polarization maps is crucial to verify the existence of 
gravitational wave background produced at the inflationary epoch. Moreover, gravitational 
lensing of clusters and hot electron scattering of reionization would be able to yield both E- 
and S-modes. To study these problems a sharp decomposition of E- and S-modes from Q 
and U maps is required. 

If both Q and U maps are available over the whole sky, one can flnd the whole sky 
maps of E- and S-modes with the spherical harmonic decomposition, because the relation 
between the maps of (Q, C/) and {E, B) in the space spanned by bases of spin two harmonics 
is local (Kamionkowsky et al. 1997; Zaldarriaga & Seljak 1997). However, the observed 
maps cannot be global; it is always limited by the contamination of our galaxy and other 
foreground sources. The relation between the maps of (Q, U) and {E, B) in physical space 
contains the Laplace operator, and therefore, it is non-local. The E/B decomposition with 
the spatially-limited maps of Q and U will not be unique if we lack of information of the 
polarization and its derivative on the boundary of the maps. 

For noiseless samples, the problem of uniqueness would be solved by constructing or- 
thogonal modes with window functions to fit the requirements of boundary conditions (Lewis 
et al. 2002; Bunn et al. 2003; Smith 2006; Smith & Zaldarriaga 2007). It is, however, similar 
to the domain (or windowed) Fourier analysis. The result will not be useful to study the 
structures in physical space (e.g. Chiueh & Ma 2002). 

The other challenge caused by the derivative operator is because the Q and U maps are 
discrete. Mathematically, the derivative operators or are continuous linear operators 
mapping functions defined in Hilbert space, while the observed samples Q and U actually 
are defined in space spanned by base Vi, i & I, which is a set of finite indices. This difference 
leads to large numerical errors when the discrete maps are noisy. 

The last, but not least, problem is from the smallness of i?-modes. On the scale of one 
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degree order, the power of B-mode caused by gravitational waves at inflationary epoch is less 
than that of £^-modes by a factor of at least 10^. As the maps of Q and U are random fields, 
the variance of the random fields will lead to the mixing of E- and S-modcs. Consequently, 
E- and B-modes would be contaminated from each other. Therefore, it is difficult to recover 
the power of i?-mode if the powers of E modes are much larger than that of B modes. 

In this paper, we develop an algorithm of the E/B decomposition based on the discrete 
wavelet transform (DWT) analysis, which is a compromise between the decompositions in 
physical space and scale space. The DWT analysis of the CMBR temperature fluctuation 
maps has attracted much attention in the last decade (Pando ct al. 1998; Sanz ct al. 1999; 
Mukherjee et al. 2000). Besides these points, we especially take the advantage of the so-called 
wavelet-Galerkin discretization (e.g. Louis et al. 1997), which is to approximate derivative 
operator to a matrix in space spanned by wavelet bases. For some available wavelets, the 
matrixes are exactly diagonal in scale-space, and narrowly banded in spatial space. This 
made the uncertainties from boundary, noises and variances are under control. We will 
study the conditions, under which the information of small B-mode can approximately be 
extracted from noisy maps of Q and U. 

The paper is organized as follows. Section 2 presents the method of the E/B sepa- 
ration in the DWT space. Section 3 tests the DWT algorithm with samples with known 
spatial structures. We show that the method effectively to suppresses the uncertainties from 
boundary effect and noise. It is also effective to identify spatial structures of E and B 
flclds. Section 4 presents the tests for samples of Gaussian random fleld. The effect of the 
variance of Gaussian random field is analyzed, especially the problem of the mixing of E- 
and S-modes. Section 5 addresses the effectiveness of the wavelet-Galerkin discretization. 
Finally conclusions are given in Section 6. The DWT representation of derivative operators 
are given in the Appendix. 



Let us consider polarization samples in a patch of sky, which can be approximated as a 
plane described by Cartesian coordinates (x, y). In this case, the fields of E{x, y) and B(x, y) 
are related to the maps of Stokes parameters Q{x,y) and U{x,y) by (Seljak 1997) 



2. 



Method 



2.1. E/B separation in DWT representation 



V^E{x,y) 
V'B{x,y) 



idl-dl)Q{x,y) + 2d,dyUix,y) 
2d,dyQ{x,y)-{dl-dl)U{x,y) 



(1) 
(2) 



-4- 



where is 2-D Laplace d^ + dy. 

We first take a wavelet-Galerkin discretization of equations (1) and (2) to rewrite these 
equations in the DWT space. We can assume that the patch is a L x L square. The size of 
each pixel is J being a integral, one can project the maps into the DWT space by 

^h,h^ J Q{^^y)(pj,h{^)(pj,h{y)d^dy (3) 

^Yi,i2 ^ J ^(^' y)<l>J,h {x)<f>J,h {y)dxdy (4) 

where (pj^ilx) is orthogonal scaling function on scale J (e.g. Fang & Thews 1998). It is 
non-zero mainly in the cell in x-space from to (/ + l)L/2'^ . The index / runs from 

to 2"^ — 1. It spans the spatial range from to L. The variables e^^ i^ and ef^ actually are 
the maps of Q and U on scale J. Since the observed maps of Q and U are always pixelized, 
the projection of eqs.(3) and (4) does not lose information if the size L/2"^ is the same as 
that of pixels of observed samples. 

One can further take a projection on eqs.(l) and (2) as 

{V^E{x,y),ct^j,i,{x)ct^j,iM) = {{dl-dl)Q{x,y) + 2d,dyU{x,y)AjM{^)My)) (5) 
{V-'B{x,y)Aj,w{^)ct>j,iM) = {2d,dyQ{x,y) - {dl- dl)U{x,y),cl>j,i,{x)ci>j,iM) ■ (6) 

With the DWT decomposition of E{x,y) and B{x,y) 

E{x,y) = ^ef;,i^<Pj,hi^)'Pj,i2{y), ^h,i2 = J E{x,y)(pj^i,{x)<j)j^i^{y)dxdy (7) 

B{x,y)^^ef^i^^j^i^{x)^j^i^{y), ef^i^^ / B{x,y)(t)j^i^{x)(t)j,i^{y)dxdy (8) 
eqs. (5) and (6) yield matrix equations 

where the matrix Ti^^ is given by 

t[J^ = J ct>j,i{x)d:ct>j,,{x)dx. (11) 
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Obviously we can do the projection of eqs.(5) and (6) using any bases in 2-D space 
L X L. However, for proper wavelet scaling functions, the integral J (l)j^i{x)d^(l)j'^i'{x)dx are 
zero for J ^ J'. This point is important for a wavelet- Galerkin discretization (see discussion 
in §4). In this case, all quantities of eqs.(9) and (10) are on scale J, and eq.(ll) gives gives 

Ti? - (12) 

where h = 1/2^ . rfi^i, is non-zero only in a narrow band \l — l'\ < M, where M is an 
integral, depending on wavelet. For Daubechies 6 wavelet, the non-zero coefficients r"_;, are 
|/ — /'I < 4, or M = 4. The values of r/_;, and rf_i, of Daubechies 6 wavelet are listed in 
Table 1 of Appendix. 

Thus, eqs.(l) and (2) defined in continuous space {x,y) are reduced to eqs.(9) and (10) 
defined in a space spanned by orthogonal bases (t)j,h{x)(j),j^i^{y). The discretized Eqs.(9) 
and (10) are an approximation of eqs.(l) and (2). Equations(9) and (10) do not contain 
information on scales less than L/2"^. However, this discretization is reasonable in the sense 
that it does not introduce false correlations, or lose information of discrete datasets Q and 
U . The derivative operator on a function defined in a space spanned by bases (l)j^i^{x)(t)j^i^{y) 
will yield a function in the same space. This is required by a wavelet-Galerkin discretization 
(§4). It ensures no signal to be produced on scales less than L/2'^. The eqs.(9) and (10) give 
a. E/B decomposition from observed maps Q and U. 

It should be pointed out that not all wavelets yield J-diagonal matrices hke eq.(12). 
For instance, the popular wavelet Daubechies 4 is not suitable for this discretization, as 
the matrix of derivative operator in space spanned by Daubechies 4 scaling functions is not 
J-diagonal. More discussion on the wavelet-Galerkin discretization will be given in §4. We 
will first study how to develop the algorithm of the E/B decomposition with eqs.(9) and 
(10). 



2.2. Maps of E,^,,, and Mi.^i, 
Equations(9) and (10) can be rewritten as follows 

E%i.'^);«'^2)<^^=E^i,'- (13) 

''li''2 

5] M(,,,,,);(,,,,,)e^,,, (14) 

where the matrix M is 

M(,,,,);(,. = T^^jA,,, + 5,,,. . (15) 
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If we use Daubechies 6, '^i^,i2 and B^j^^j s-^® given by 

4 4 4 

mi=—4: m2=— 4 mi,m2=— 4 

4 4 4 

■^h,«2 — / J mi-^ m2^h+mi,l2+m2 / j mi Si+mi,;2 ^ / j Hi ,l2+Tn2 ■ ' J 

mi,m2=— 4 mi=— 4 m2=— 4 

The equations (13) and (14) look like the matrix equations of the DWT variables of E 
and B fields, ef^ i^ and e^^ ^^- Equations(16) and (17) give, the sources '^1^,12 and 'Bi^,i2 on the 
right hand side of eqs.(13) and (14), respectively. It seems that one can separate E/B by 
solving the matrix eqs.(13) and (14). However, with the coefficients t/^^ given in Appendix, 
we can show 

J2^ihM,(i[,i'2) = ^- (18) 

That is, the matrix M is singular. One cannot use a standard linear solver to solve cqs.(13) 
and (14). This problem, of course, is directly related to the non- uniqueness of the solutions 
E{x, y) and B{x, y) given by the Poisson equations (1) and (2) without knowledge of bound- 
ary conditions. We will not try to solve the singular matrix equations (13) and (14), but 
directly use eqs.(16) and (17) for the E/B decomposition. 



2.3. E/B decomposition with Ej^^jj and ^1^,12 

The spatial resolution of the source terms E^^^^j and B^^^^j is the same as maps Q and 
U. It can be used to calculate the DWT power spectrum of V'^E and V^S fields, and other 
statistics. To do these, we should first find the wavelet function coefficient (WFC) of the 
maps ¥.i^^i2 and B^^^^j by 

?i = ECj,M'Ei' (19) 
1' 

€1 = E ^J'M'Bi' (20) 
1' 

where, for simplification, we use 2-D vector notation defined by j = (71,72) and 1 = {hjh), 
and li = 0...2^^ — 1, I2 = 0...2^^ — 1; ji, 72 can be any integral less than J. Index (j, I) refers 
to the cell on scale j and at position 1. The 1x1' matrix Cj,i,i' is given by 



Cj;(i,i/) = / i^juhi^)i^j2,i2{y)(pj,i{{x)(pj,i'2iy)^^^y 



(21) 
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where ipj,i{x) is 1-D wavelet function referring to cell on scale j and at position I. Cj^(i^i/) is 
a banded matrix with respect to 1, 1'. Therefore, the relation between e?i, e?i and Q, U are 
spatially quasi-local. 

With the WFCs, the DWT power spectrum is given by (Fang & Feng 2000) 

^E,B ^ ^^^E,B^2^ ^22) 

where (...) is the average over all cells 1. One can directly use the DWT power spectrum to 
measure E- and S-modes. Pj ' is banded Fourier power spectrum. For a statistically homo- 
geneous random field, the DWT power spectrum is related to the Fourier power spectrum 
P^'"(ni,n2) of (E, B) maps by 

^ oo 

p_E,B^^ E |V'(ni/2^-^)V'(n2/2^-^)|2p«'»(ni,n2). (23) 

ni,n2=—oo 

Clearly, Pj is banded Fourier power spectrum with the window function 

Wi{ni,n2) = -^|V^(ni/2^'i)V^(n2/2^'^)|l (24) 

Function ip(n) is the Fourier transform of the basic wavelet. Pj contain all valuable quantities 
of second order statistics from random samples in a finite area L x L and pixel . The 
window function may cause spurious features and false correlation, such as aliasing effect, in 
the Fourier power spectrum. With the DWT analysis, the ahasing effects can be effectively 
suppressed (Fang & Feng 2000). 



2.4. Effect of noise on power spectrum 

The maps of Q and U are usually noisy and can be given by Q + AQ, and U + AU. 
The DWT variables of noisy DWT maps are then + AQi and eY + AUi, where AQi 
and AU\ are the DWT variables of AQ and AU. Assuming the noise is Gaussian and 
statistically homogeneous, the DWT variables AQi and AUi of noise have to satisfy the 
following statistical properties 

{AQAQi') = ct'q6i,v, {AUAUv) = a'uS,,r, (AQAUy) = (25) 

where ag and au are the variance of the noise of Q and U maps, respectively, and are 
independent of I. 
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Using AQi and AC/i to replace and in eqs.(16) and (17), we can construct the 
noise maps of AEi and ABi. First, with eq.(25) we can show 

(AEiABi,) = 0. (26) 

That is, noise eq.(25) does not cause false correlation between E and B modes. This is 
very helpful for the E/B decomposition. Second, with eqs.(16) and (17), one can find the 
variance of the noise maps AEi and ABi to be 

((AE,)') = N^^^afj + iV^'Vj (27) 

((ABi)2) = Ar(i)(7j + Ar(2)a2 (28) 

where 



m=l 



(1)^2 
m J 



2)12 



(29) 



m=l 



N'^^^ and N'^'^^ are from the terms containing T^^], and T^%, respectively, in eqs.(16) and (17). 
For Daubechies 6 wavelet, VNW ~ 1.2 and Vn^ ~ 5. That is, in the Daubechies 6 DWT 
algorithm [eqs(16) and (17)], the operator of derivative d does not significantly change the 
level of the noise, while the operator for leads to an increase of the variance by a factor of 
5 with respect to the variance of AQ and AU maps. This shows that derivative will generally 
amplify the effect of noise. However, the matrix t/J^^ is exactly diagonal with respect to j, 
the derivative operator in the DWT representation does not transfer the noise from one scale 
to others. In this sense, we have a handle on the noise. 

As noise and signal are statistically uncorrclatcd, the power spectrum of Ei and Bi can 
be reconstructed by subtracting the power of noise as 



pE ^ pE 
j j 

pB _ pM 
j j 



Pi 



AE 



(30) 
(31) 

where and P^ are the DWT power spectrum of maps Ei and Bi given by eqs.(16) and 
(17), respectively, using noisy Q and U. P-^^ and P-^^ are the DWT power spectra of noise 
maps AEi and ABi. The algorithm of subtracting the noise DWT power spectrum P.^^ and 
P^^ scale-by-scale is similar to the subtraction of shot noise power from the DWT power 
spectrum of galaxy survey (Fang & Feng 2000). 



2.5. The effect of boundary 



For a sample of finite area, the DWT power spectrum analysis does not need a window 
function to treat the spatial domain. The effect of boundary can effectively be reduced by 
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dropping the DWT variables related to cells {j, I) located on or near the boundary (Pando & 
Fang 1998). When derivative operators, d^, dy, are involved, the boundary effect would be 
more serious, because the matrix of derivative operators in the DWT representation is not 
exactly diagonal with respect to the spatial index /. Nevertheless, the matrix T^^} [eq.(ll)] is 
narrowly banded, the effect of boundary can still be reduced by dropping boundary modes. 



To test the DWT algorithm developed in §3, we consider, in this section, samples with 
given spatial structures, and compare the maps E; and given by cqs.(16) and (17) with 
that directly calculated from E and B. This comparison is only to test the discretization of 
derivative operator, but says nothing about the amount of information loss associated with 
the algorithm. 



We use two scalar functions ipE{x,y) and ipB{x,y) to produce E and B maps in 2-D 
space by the following way 



One can then produce the DWT variables ef^; and efi, with eqs.(7) and (8). With these 
results, we can further produce the maps of Ki^i/ and B;^// with eqs.(13) and (14). Thus, for 
given scalar functions iPe{x, y) and ipsi^, y), we have the samples and B; and then, the 
the DWT power spectrum and other statistical properties of maps E; and B; The ratio 
between the powers of E- and S-modes can be adjusted by the ratio between the functions 
■0B and i/jb- 

On the other hand, using the function ipE^x^y) and ipBi^jy), one can produce the 
samples of the Stokes parameters Q{x,y) and U{x,y) maps by 



3. Tests with samples having known spatial structures 



3.1. Samples 



E = -V'V^E, B = -VVs- 



(32) 



Q{x,y) 

U{x,y) 



{dA - dydy)ijjE{x, y) - 2d^dyijjB{x, y), 
2d^dy'il)E{x, y) + {dA - dydy)'il)B{x, y). 



(33) 
(34) 



We add Gaussian white noise in the Q and U maps pixel-by-pixel with signal-to-noise ratio 
equal to 10, 50, and 100. These Q and U maps are used as the simulation of observed 
samples. 
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With noisy maps Q and C/, we can produce the variables ^, and e[^, by the projection 
of eqs.(3) and (4). Finally, we have maps E;^;/ and B;^;/ using eqs.(16) and (17). Thus, wc 
can test the algorithm by comparing the statistics of the maps ;/ and B; given by Q and 
U [eqs.(33) and (34)] with that directly calculated from E and B of eq.(32). 



3.2. Recovery of spatial structures 

The scalar functions ipE and -05 are taken to be sample A.) Gaussian function ipE,Bix, y) — 
aE,B —{x"^ +l/^)/2'i^; sample B.) the Lcgcndrc function ipE,B{,x,y) = aE,BPm{x)Pm{y). 
Both samples are in the area —0.5 < a; < 0.5 and —0.5 < y < 0.5 and pixels 512 x 512, i.e. 
J = 8. The coefficients aE and are used to adjust the ratio of the powers of E^ and 
B/^//. Wc use aE = I and as = 1/10. That is, the power of ii^-mode is larger than B mode 
by a factor 10^. The maps of E;^;/ for samples A and B in the central square 32 x 32 pixels 
are shown Figure 1. The maps of B; have the same shape of Figure 1, but the intensity is 
weaker than Figure 1 by a factor 10. 




Fig. 1.— The DWT maps of Ei of sample A: ipE = cxp[-{x^ + y'^)/2d\ and 2d^ = 1600 
(left), and sample B: -0^ = iI^b/clb — Pm{x)Pm{y), and m — 100 (right). 



As mentioned in §3.1, with ipE and ips one can produce the maps of Q and U by eqs.(33) 
and (34). Using noise-added maps of Q and U, we can further calculate noisy variables E^^^^j 
and B^^^jj with eqs.(16) and (17). This is the recovered maps of E^^^jj '^hh- Figures 2 
and 3 present the recovered maps of E;^^;^ and B^^^^j for samples A and B, respectively. The 
signal-to-noise ratios (S/N) are S/N=10, 50 and 100 from left to right. Comparing Figures 2 
and 3 with Figure 1, we can conclude that the spatial structures of E mode maps can be well 
recovered with the noisy maps Q and U if S/N > 10. The recovery of B mode structures 
is relatively poor. For sample A (Fig. 2), we may pick up the original structures of B field 
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Fig. 2. — The maps IEzi,/2 (^P panels) and B;^ ^2 (bottom panel) of sample A, which are 
recovered by eqs.(16) and (17) with noisy maps Q and U, and the signal-to-noise ratios are 
taken to be S/N=10 (left), 50 (middle) and 100 (right). 

with all S/N> 10 noisy maps of Q and U, while for sample B, the structures of B field can 
not be seen with S/N — 10 map. That is, the structure identification of sample A is much 
better than sample B. This is because the field of sample A is highly inhomogeneous, while 
sample B is not so inhomogeneous. The latter is easily contaminated with a statistically 
homogeneous Gaussian field. 

3.3. Recovery of power spectrum 

We now turn to the recovery of the power spectrum. For wavelet Daubechies 6, the 
non-zero elements of the matrix T^, are in a band \l — l'\ < 4, and therefore, cells distant 
from boundary, larger than Al — 4, will be less affected by the boundary. Thus, one may 
expect that the power spectrum recovery would be reasonable with dropping 4 boundary 
cells. 

Figures 4 and 5 show the DWT power spectra of E;^;/ and B^^^/ of both original and 
recovered samples of sets A and B, respectively. It includes 1.) the power spectra of the 
original maps, i.e. the map directly given by E and B [eq.(32)]; 2.) the ratio P^/PJ, where 
P^ is the original power spectra from maps E^^^^j and B;^^;^ of eq.(32), and P^ is the recovered 
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Fig. 3. — The same as Figure 3, but for sample B. 

power spectra of maps Q and U with eqs.(16) and (17) without dropping boundary cells; 3.) 
the same as 2.) but dropping 4 boundary cells in the recovered power spectra. 

In Figures 4 and 5, the scale j = (ji, 72) is described by an effective scale defined as jcff 
as (l/2-^<=*f) = [(1/2-^^)^ + (l/2-'2)^]^/^. The samples arc symmetric with respect to x ^ y. The 
power of mode (ji, should be the same as {j2,ji)- Thus, for J = 8, the available pairs 
(j2,Ji) are (4,4), (4,5), (4,6), (4,7), (5,5), (5,6), (5,7), (6,6), (6,7) and (7,7), corresponding 
to Jeff = 3.50, 3.84, 3.95, 3.98, 4.50, 4.84, 4.96, 5.50, 5.84, 6.50. The modes on scales with 
ii, J2 < 3 are dropped, as all cells are affected by the boundary effect. 

We see from Figures 4 and 5 that the power spectra can indeed be well recovered by 
dropping 4 boundary cells. However, the boundary effect of sample A are less serious than 
sample B. This is because, for sample A, both iPe,b{x, y) and dnipE,B{x, y) are very small at 
boundary. The contribution to power by boundary cells is low. On the other hand, for sample 
B, the power spectra recovered without dropping boundary cells are significantly different 
from the original one. The i?-mode power spectrum is hugely affected by the boundary. On 
small scale the derivative operator in the DWT representation is determined by data at a 
few discrete points, which leads to large error. This problem is always present in algorithms 
involving taking derivative on discrete data sets. Nevertheless the error caused by boundary 
is decreases rapidly as the scale increases. 

To measure how good the recovery of power spectrum is, we use the ratio P^/P[ to 



Fig. 4. — Left panel: The DWT power spectra Pj of the original maps Ei and Bi of sample A. 
The power spectrum on the top is for the E mode, and the lower is for the B mode. The ten 
data points correspond to (ji,j2) = (4,4), (4,5), (4,6), (4,7), (5,5), (5,6), (5,7), (6,6), (6,7) 
and (7,7) from left to right. The power ratio E/B is equal to 10^. Middle panel: the ratio 
P?/PT, where P? is the original power spectra from maps E^^^^j B^^^^j of eq.(32), and 
PT is the recovered power spectra of Q and U with eqs.(16) and (17) and without dropping 
boundary cells. The top is for E mode, and bottom is for B mode. Right panel: the same 
as middle panel, but with 4 boundary cells dropped. 

describe the deviation of the recovered power spectrum with noisy maps from the original 
one. All error bars are the variances calculated from 100 independent noisy maps. The 
results are plotted in Figures 6 and 7. 

First we see that the effect of Gaussian noise is small at larger scales, because the noise 
is added on each pixels (finest scale) of the maps Q and t/, and the uncertainty on large 
scales is suppressed. This point can also be seen from the fact that the error bars of modes 
(4,7), (5, 7), (6, 7) and (7,7) are much larger than others. It is because the Gaussian noise on 
smallest scales, ji or j2 — 7, is not suppressed. 

Figures 6 and 7 show that other than the modes with ji or j2 = 7, the power of E mode 
can be reasonably recovered up to mode (6,6), or Jeff = 5.5, when S/N= 10. As expected, 
the recovery for B mode generally is poor. Nevertheless, we can recover the DWT powers of 
B mode till (5,5), or Jeff = 4.5, when S/N=50 (sample A) or S/N=10 (sample B). 

An interesting point shown in Figures 6 and 7 is that the effects of noise on samples A 
and B are different. The error bars of sample A generally are larger than that of sample B. 
This is probably because for sample A, other than the central part, most cells are smooth, and 
have low local fluctuations. For those cells, the fluctuations of noise will strongly contaminate 
the power of original fleld, especially when derivative is involved. On the other hand, for 
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Fig. 5. — The same as Figure 4, but for sample B. 

sample B, most cells have relatively stronger local fluctuations, and the effect of noise is 
relatively low. 

4. Tests with samples of Gaussian random fields 
4.1. Samples 

With the preparation given in the previous section, we can consider the case that 
ipE{x,y) and ipB{x,y) as random flelds. The sample of E and B can still be generated 
with the same procedure of §3.1, but il)E{x,y) and i/jB{x,y) are taken to be Gaussian ran- 
dom fields with Fourier power spectra aEk'"^ and ask""", respectively, and a = 3.6. The 
variable k"^ = + k^, and k^, ky arc the Fourier variables of y space, respectively. The 
constant factors and as are used to adjust the ratio oi E/B power. From eq.(32), the 
power spectra of E and B are 

Pe = aEk''-", Pb = aBk''-". (35) 

We produce the maps in an area described by coordinate {x, y) in range —0.5 < a: < 0.5 and 
—0.5 < y < 0.5 with pixelized into 512x512. With maps E and B, one can find the maps E 
and B by eqs.(13), and (14). 

The DWT power spectrum of E and B is shown in Figure 8. Since J = 8, the available 
modes (ji, j2) still are (4,4), (4,5), (4,6), (4,7), (5,5), (5,6), (5,7), (6,6), (6,7) and (7,7). The 
modes with ji, j2 < 3 are dropped, as they have only boundary cells. The ratio of the E/B 
power is taken to be 10. The error bars are from the variance of 100 samples. We see from 
Figure 8 that the powers of modes (4,7), (5,7) and (6,7) are nearly about the same as (7,7). 
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Fig. 6. — The ratio Pj/PJ for sample A, where P" is the power spectrum of the original 
maps IEii,i2 ^^1,^2 5 Pj is calculated by eqs.(30) and (31) from Q and U with noise 
addition on the level of S/N equal to 10 (left), 50(middle), and lOO(right). In each panel,the 
top is for E mode, and the bottom is for B mode. 



Similarly, the powers (4,6), (5,6) are nearly about the same as (6,6); the power (4,5) is nearly 
about the same as (5,5). It is because in the case of ji < j2, the power is dominated by the 
small scale j2. 



4.2. Effects of random field 

Unlike the maps in §3, all the maps of ipEix,y), ipBix,y); Q, U; and E, B are random 
fields. A serious problem caused by random fields is that the power of E mode may leak to 
B mode, and vice versa. That is, even when the original B mode power is zero, the recovered 
B mode power would not be zero. 

To demostrate the power leakage, we take iPe{x, y) to be a Gaussian random field with 
the same Fourier power spectrum as Figure 8, while ipB{x,y) to be 0. That is, the power 
of B mode originally is zero. Figure 9 presents the recovered DWT power spectra of E- 
and i?-modcs. We see that the recovered £'-mode power spectrum is nearly the same as the 
original one shown in Figure 8, except that the recovered power spectrum on the finest scale 
is a little smaller than the original one. However, the recovered i?-mode power spectrum is 
not zero. It is spurious B power. It arises from the leaking of E mode power to B mode. On 
large scales jes < 4, the ratio of the E/B power is about 10^, while on small scales jes > 6, 
this ratio is less than 10^. This is caused by the variance of random field. Thus, one may 
conclude that for the 512 x 512 sample of a Gaussian random fields of eq.(35), the developed 
algorithm would be effective only if the ratio E/B is less than 10^ on small scales. 
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Fig. 7. — The same as figure 6, but for sample B. 



As a comparison, we plot Figure 10, in which the ipE{x,y) is given by samples A and 
B, while y) = 0, i.e. the power of B mode originally is also zero. Figure 10 shows 

that the recovered B-mode powers are also not zero. However, it generally is less than the 
original one by at least 3 orders. That powers seem to come from the numerical processes. 
Therefore, the errors caused by the variance of random fields are serious. 



4.3. Recovery of E, B power spectra 

As in §3.3, we measure the soundness of the recovery of power spectrum by the ratio 
Jj^/Jj", where P-* is the power spectra of original maps E;^^;^ and B^^^^j from eq.(32), and P-' 
is the recovered power spectra from noisy maps of Q and U. The Gaussian noise added on 
the maps Q and U are on the levels S/N=100, 20 and 10. Similar to §3.3, four boundary 
cells are dropped. The results are plotted in Figures 11, 12 and 13, for which the ratio of 
the powers of E and B are equal to, 10, 20 and 100, respectively. 

Figure 11 shows the powers of E mode can be recovered on all scales jeff < 6 on all noise 
levels. On the smallest scale, jcs = 6.5, the ratio Pj/Pj of E mode is slightly lower than 1. 
This deviation is almost independent of the level of S/N. Therefore, the errors mostly are not 
due to the Gaussian noise addition, but from the effect of leakage. This point is consistent 
with the leaking shown in Figure 9, which also give a little small power on the smallest scale. 
More interesting. Figure 11 shows that the B mode can be perfectly recovered on all scales 
and all noise levels considered. 

Figure 12 presents the case of E/B = 20. The results on scales jcff < 5.5 are about the 
same as the case of E/B = 10, while the deviations of from Pj on scales jefr > 5.5 are 
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Fig. 8. — The power spectra of E (higher one), B (lower one), for which ipE and ipB are 
Gaussian random field with the Fourier power spectra aEk~^'^ and ask'^'^, respectively. 
The ten data points correspond to (ji, j2) = (4,4), (4,5), (4,6), (4,7), (5,5), (5,6), (5,7), (6,6), 
(6,7) and (7,7) from left to right. The error bars are from 100 samples of the Gaussian 
random field ipEi^, y) and iPb{x, y). The ratio of powers of E and d B modes is equal to 10. 

larger than that ol E/B — On small scales, the recovered E powers, PJ, are little lower 
than the original power P?, while the recovered B powers are little higher than the original 
power. 

Figure 13 is for the case of E/B — 100. It shows that the recovered E mode power 
spectrum is still good on scales jeff < 5.5 for all S/N. However, on scales j^s > 5, the 
recovered B mode power spectrum generally is higher than the original one. This deviation 
is expected, as Figure 9 shows that the leaked power from E mode to B mode can be as high 
as 1% on small scales. Nevertheless, the recovered B mode power spectrum is reasonable on 
scales Jeff < 5, even when the S/N is equal to 10. That is, one can pick up the weak signal 
of B mode with the DWT algorithm even when the Gaussian noise level of Q and U maps 
is comparable with the B mode signal. 



5. Discussions 

The relationships between the polarization maps of {E, B) and (Q, U) are differential. 
For pixelized samples of Q and U in finite area, the algorithm oi E/B decomposition should 
be able to properly handle the derivative operation on a spatially discrete and noisy data 
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Fig. 9. — DWT power spectra of E mode (top) and B mode (bottom). The iljE{x,y) is a 
Gaussian fields witli tlie same Fourier power spectrum as Figure 8, wliile taking ipsi^, u) — 0. 

set. The derivative operator is a continuous hnear operator to mapping functions defined 
in Hilbert space, while the functions Q and U are defined in space spanned by bases Ui, in 
which i is a set of finite index. Therefore, we should approximate the derivative operator 
from mapping between functions defined in Hilbert space, to a mapping in subspace spanned 
by Hi. 

What we need to calculate is 



where O is a linear continuous operator, like Laplace or derivative, and / and g are function 
of X, y in continuous space < x,y < L. However, we don't know /, but only the discretized 
/, which is given in N pixels (cells). That is, / can be expressed as 



k=l 

where function Wi{x,y) is the binning function of pixel k, and ak is the observed / at pixel 
k. The simplest binning function would be the top-hat window function of pixel k. 

With A^-dimensional space spanned by bases [vi{x, y), ...vn{x, y)], eqs.(36) and (37) yield 



Of = g, 



(36) 



N 




(37) 



N 

9i = {9,Vi) = ^{Owk,Vi)ak. 



(38) 
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Fig. 10. — DWT power spectra of E mode (top) and B mode (bottom). The ipEix^y) is 
given by sample A and sample B, while taking iJjb{x, y) = 0. The left panel is for sample A, 
while the right panel is for sample B. 

The matrix O^k = {Owk, Vi) gives a discretization of operator O from the space Q < x,y < L 
to a finite-dimensional subspace [f . 

A Galerkin discretization requires the following equation to be hold for all vi 

{{g-Of),Vi)^0 (39) 

That is, eq.(36) should be hold in the subspace spanned by bases [vi{x,y), ...VN{x,y)]. In 
this case, the matrix Oik is a linear operator to map functions defined in the subspace Vi. If 
/ is a function of the subspace Vi, g = Of is also a function of the subspace. 

It can be seen from eqs.(36) - (38) that the discretization of operator O actually is 
inevitable for all algorithms. To treat the data eq.(37), we must use some base Vi in the 
spatial domain. It will yield a matrix Oik, regardless whether eq.(39) is hold with the bases 
Vi- The Galerkin method gives a best discretization of eq.(36) or the operator O (e.g. Louis 
et al. 1997). 

We use (})j_i{x)(t)jji{y), 1,1' = 0...2"' — 1 to be the bases to span the subspace with 
dimension = 2"' x 2*^. It can be shown that the conditions of Galerkin discretization, 
eq.(39), will be satisfied for operator O — d^, dy, and dy. That is, for any function f{x, y) 
of the subspace spanned by bases (t)j^i{x)(t)j^i'{y), 1,1' — 0...2"^ — 1, the result of Of are also 
functions of the subspace. This is the wavelet-Galerkin discretization. Matrix {Owk,Vi) will 
be invertible. In this sense, the discretization does not lose information, or introduces false 



Fig. 11. — The ratio PJ /Pj for Gaussian random field iPe,b{x, y) with Fourier power spectra 
aE,Bk~^'^. The top panels are for E mode, and bottom for B mode. The ratio of the powers 
E/B is equal to 10. The maps Q and U are added noises with the level of S/N equal to 100 
(left), 20 (middle) and 10 (right). 

data or correlations. 

Obviously, the Galerkin discretization is not unique. One can use different wavelets to do 
the Galerkin discretization. To apply the discretization, the matrix Oik = {Owk, f i) should 
have the following desirable properties. First, the matrix Oik has to be sparse, narrowly 
banded. In this case, one can effectively minimize the information lose due to dropping 
boundary modes. A narrowly banded matrix can also effectively reduce the spreading of 
errors among cells with different 1. Secondly, in order that the errors not increase with the 
size of the matrix (number of data), the "width" of the band in which the matrix elements 
Oik is non-zero, should be independent on N. 

6. Conclusions 

The algorithm developed in this paper can be summarized as follows 

• From observed noisy and discrete maps Q{x,y) and U{x,y) we calculate their DWT 
maps Qii,i2 and Ui-^,i2 on the finest scale j = {J, J), which is given by the resolution. 

• Using eqs.(16) and (17) we decompose Qii,i2 and Ui^^j^ i^^o ^i^^i^ and B^j^^^. 

• Using eqs.(19) and (20) we calculate WFCs e^'j and e^j on scales (^1,^2) and ji,j2 < J- 

• Using the WFC maps if^ and if^ we calculate the DWT power spectra by eqs.(30) and 



Fig. 12. — The same as Figure 11, but the power ratio of E and B modes is equal to 20. 



(31). 

• We identify spatial structures with maps of ^i^j,^ and Bjj 

With this algorithm, it is possible to recover the power spectrum of S-mode random 
fields from noisy Stokes parameter maps Q and U when the power ratio E/B is as high 

as 10^, and the S/N is equal to or higher than 10. For samples with given structures, the 
i?-mode structure can also be identified when the power ratio E/B is equal to 10^. Besides 
power spectrum, the DWT variables of SFCs (e^'j, e^j) and WFCs (e^'j, e^j) can be used for 
high order statistics, such as high order moments, scale-scale correlation, cross correlation 
between the E and B and other maps. 

With the DWTs, one can construct orthogonal, divergence-free vector wavelets. It has 
been used for a local analysis of the velocity field of incompressible turbulence (Urban 1995; 
Kishida et al. 1999; Albukrek et al. 2002). The divergence- free B field is similar to a 2-D 
velocity field of turbulence of incompressible fiuid (e.g. Pina 1998). Therefore, it would 
be valuable to further study the DWT E/B decomposition with the divergence- free vector 
wavelets. 
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Fig. 13. — The same as Figure 11, but the power ratio of E and B modes is equal to 100. 



A. Derivative operator in wavelet representation 

In the DWT space, the operators of derivatives are represented as a matrix 

t\Pj,3' = J hi(^)d:<l>rA^)dx. (Al) 

That is, the matrix is diagonal with respect to rjj], is given by (Beylkin 1992; Kwon 
1998) 

- l^n (A2) 

where h = 1/2K The matrix elements r["|, depend on the type of wavelet. For Daubcchies 6 
wavelet, the non-zero coefficients r|'"], arc \l — l'\ < 4. The values of rj-"], are listed in Table 
1, in which m = I — I' . Therefore the coefficients of Tm^ of eqs.(ll) and (12) are given by 

Tln^ = • (A3) 

It is interesting to see that the non-zero band of first and second order derivative oper- 
ators dx and are the same. This is different from the estimation of derivative operator by 
differential approximation. 
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